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Introduction 


The internal flow field plays an important role in numerous aspects of the 
design and operation of a solid propellant rocket booster, such as is used in 
the current shuttle. During design, accurate predictions of the internal flow 
field are needed to properly estimate motor case insulator requirements, and 
thereby avoid performance penalties due to excessive inert weight. Under- 
estimation of insulator thickness has an even greater penalty, as it can result 
in burnthrough of the motor case and failure of the booster during flight. 
Pathological flow problems in slots, seals and around inhibitors can create 
erosion and ablation problems which must be eliminated, and internal flow 
predictions can give valuable insights into these potential problems. 

The capability for predicting the internal flow field and the associated 
particle history in the rocket motor is also essential to motor performance, 
nozzle erosion, motor stability and modeling of the erosive burning 
characteristics of solid rocket propellants. Analysis of the rocket's internal 
flow field using the techniques of Computational Fluid Dynamics (CFD) is 
frequently the only feasible means for evaluating other motor phenomena not 
measurable through static motor tests. 

Predictive techniques used commonly at present involve the use of inviscid 
methods which will be obviously inadequate as the flow transitions to turbulent 
regime. Techniques estimating the chamber mean flow from experimental data and 
correlations of one- or two-dimensional inviscid core flow calculations with 
boundary layer corrections to the mean flow estimated iteratively require the 
division of the flow into an inviscid region and a boundary layer. In most 
cases of practical interest, such as in the vicinity of the slots at the 
segment joint, such a division is not possible and the boundary layer 
approximations are not valid. Hence, a more general approach for flow field 
prediction is needed. 

Under Contract F04611-83-C-0003 with the Air Force Rocket Propulsion 
Laboratory, a program was initiated to accurately model the core and near wall 
regions of the rocket internal flow field, including slots, inhibitors and the 
nozzle. An existing SRA Navier-Stokes computer program (MINT) has been used to 
simulate the two-dimensional (axisymmetric) steady state flow of a nonreacting 
ideal gas in a cylindrical port cold flow model (Refs. 1-2). The computer 
program was modified to model the massive wall-blowing necessary to simulate 
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the surface of a burning solid propellant. Geometric capabilities of this 
computer program were enhanced to accomodate rocket motor geometries with 
radial slots and restrictors. Preliminary computed results show good agreement 
with the available experimental data from cold flow experiments on the axial 
velocity profile in the near wall region. This effort has demonstrated the 
value of Navier-Stokes predictive models for application to solid propellant 
rocket motor internal flows. Under the present contract, this code was 
extended to three dimensions and applied to the simulation of axisymmetric as 
well as asymmetric flow field in the vicinity of the aft field joint in the 
space shuttle solid rocket booster (SRB). 

Analysis 

The internal flow field in the solid rocket motor chamber was simulated by 
numerical solution of the ensemble-averaged Navier-Stokes equations in a 
body-fitted, cylindrical polar coordinate system. The equations to be solved 
are the continuity, momentum conservation, energy conservation and the 
turbulence model equations. The governing partial differential equations were 
formulated in conservation form by application of a Jacobian transformation to 
the equations in a cylindrical-polar coordinate. An outline of the 
transformation as well as the transformed system of equations have been 
presented in Ref. 1. The vector form of the equations is described below. 

The continuity equation is written as 

— + V.pU = 0 (1) 

at 

The momentum conservation equation is 

8(pU) + 7*(p5?) = -Vp + V*T (2) 

at 

where t is the stress tensor (molecular and turbulent) written in cylindrical- 
polar coordinates. 

The energy conservation equation is written as 
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where $ is the viscous dissipation per unit volume and q is the heat flux 
vector. The stress tensor, given by 


T ij 2 y eff e^j - — y e ff V*U 5ij 


(A) 


where the rate of strain tensor, e^j , is given by 

1 /Su.- , Su4 

e ii = — l — ~ + — 1 
^ 2 3xj 3xj 

The dissipation rate, 4>, is expressed as 


* * vi e ff [2e i je i j - (7*U) 2 ] 


( 5 ) 


( 6 ) 


The effective viscosity, P e ff, is the sum of the molecular and turbulent 
viscosities 

U e ff = V + Ut (7) • 

The turbulent viscosity, p-p* obtained from the turbulence model. The heat 
flux vector, q, is given by 

q = — (k + k t ) VT (8) 

where < and top are molecular and turbulent thermal conductivities, 
respectively. In the present analysis, k and Kf are obtained assuming 
constant molecular and turbulent Prandtl numbers Pr and Prj, i.e. 
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The turbulence model used under the present effort is a transitional 
two-equation (k-e) model discussed by Jones and Launder (Ref. 3), which is a 
modified form of the k-e turbulence model originally developed by Harlow and 
Nakayama (Ref. 4). The modifications devised by Jones and Launder allow the 
k-e model to be utilized throughout the viscous sublayer without any wall 
function assumption. Further details of the turbulence model are available in 
Refs. 1 and 3. 

The numerical procedure used to solve the governing equations is a 
consistently split linearized block implicit (LBI) scheme of Briley and 
McDonald (Ref. 5). The method centers around the use of a formal linearization 
technique adapted for the integration of initial value problems. The governing 
equations are replaced by an implicit time difference approximation. Terms 
involving nonlinearity at the implicit time level are linearized by Taylor 
series expansion in time about the solution at the known time level, and 
spatial difference approximations are introduced. The result is a system of 
multidimensional coupled linear difference equations for the dependent 
variables at the implicit time level. To solve these difference equations, the 
Douglas-Gunn (Ref. 6) procedure for generating alternating direction implicit 
(ADI) schemes as perturbations of fundamental implicit difference schemes is 
introduced. This technique leads to systems of coupled linear difference 
equations having block-banded matrix structure which can be solved efficiently 
by standard block elimination methods. 

The linearization technique permits the solution of coupled nonlinear 
equations in one space dimension by a one-step noniterative scheme. Since no 
iteration is required to compute the solution for a single time step and only 
moderate effort is required for solution of the implicit difference equations, 
the method is computationally efficient; this efficiency is retained for the 
multidimensional problems by using ADI techniques. The method is also 
economical in terms of computer storage, in its present form requiring only 
two time levels of storage for each dependent variable. Further details of the 
LBI scheme are available in Refs. 5 and 7. 
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Application and Results 


The analysis discussed above was incorporated into the MINT (Multi- 
dimensional Implicit Navier-Stokes Time -Dependent ) code and used by SRA under 
Contract F0461 1-83-C-0003 with the Air Force Rocket Propulsion Laboratory to 
simulate the steady state flow field in a cylindrical port cold-flow model. 

The code was also used subsequently to simulate the internal flow in a cold- 
flow model with a radial slot. The results of this calculation and comparison 
with experimental data are available in Ref. 2. Under the present contract, 
the MINT code was used to simulate the internal flow in the vicinity of the aft 
field joint slot of the space shuttle solid rocket booster (SRB). Calculations 
included simulation of axisymmetric flow in the chamber to estimate the 
possible over-pressure at the joint, as well as asymmetric three-dimensional 
flow in the chamber to estimate the circumferential flow and pressure 
distribution at the joint caused by pressure asymmetry at the outflow (which 
could be caused by a canted nozzle) and failure of a section of the castable 
inhibitor. 

A schematic diagram of the compuational domain for the axisymmetric case is 
presented in Fig. 1. While the region of interest is the vicinity of the slot, 
it is essential to locate the inflow and the outflow boundaries sufficiently 
far away from the slot in an attempt to keep the flow field in the region of 
interest free from the assumptions/approximations made at these boundaries. 

For this reason the inflow and outflow boundaries have been located about 12' 
away from the slot in either direction while the slot width is only 1.5". The 
calculations for the axisymmetric simulation were initiated from a quiescent 
state with pressure and temperature corresponding to the stagnation conditions 
for the chamber. The mass influx at the boundaries of the computational domain 
representing the grain surface was specified to be a product of a "permeability 
factor" and the difference between a specified "plenum pressure" and the static 
pressure at the boundary. The plenum pressure was set to equal the chamber 
stagnation pressure, thereby resulting in no mass influx at the initial 
condition. At the inflow boundary, linear extrapolation of static pressure was 
used as a boundary condition and the axial velocity at the centerline was 
computed from the chamber stagnation pressure and the static pressure at the 
centerline. The axial and radial velocities at all other locations on the 
inflow boundary were computed using the centerline axial velocity and the 


- 5 - 


velocity profiles obtained from the inviscid rotational flow analysis of Culick. 
(Ref. 8). All dependent variables except pressure were extrapolated at the 
outflow boundary. The static pressure at the outflow boundary was gradually 
lowered, resulting in mass efflux from the domain and reduction in the static 
pressure in the chamber. The reduced static pressure in the chamber started 
mass influx from the inflow boundary, as well as the boundaries representing 
the grain surface. After a flow field was established in the domain, the 
boundary condition at the grain surface was changed to specify the actual mass 
flux computed from data on grain regression rate, etc. The static pressure at 
outflow was adjusted to obtain the desired mass flow at the inflow boundary and 
the solution procedure marched in time to obtain the steady-state solution. 

For the three-dimensional calculations, the computational domain included 
the same segment length as used in the axisymmetric calculations. In the 
circumferential direction, a 180° section (-90°<^ 8 _<+90°) of the SRB chamber 
was included. The initial flow field for the three-dimensional calculations 
was obtained from the converged solution for the axisymmetric calculations. 

The appropriate boundary conditions were then changed to introduce the 
asymmetry and the solution procedure marched in time to obtain the final 
solution. 

Computational Grid 

The grid generation technique used to generate a grid of the computational 
domain is similar to that discussed by Oh (Ref. 9). Every physical location 
which must coincide with a grid point is designated as a "cluster point". 
Examples of such locations in the axial direction are the two faces of the slot 
and in radial direction they would be the grain radius, the beginning of the 
inhibitors, etc. Complementary error functions are used to distribute the grid 
points subject to the constraints of the cluster points. Figures 2 and 3 show 
the grid spacing in the radial and the axial direction, respectively, plotted 
in the computational space. Since it is necessary to resolve the flow field 
near the solid surfaces as accurately as possible, a high degree of mesh 
stretching has been used. 

These figures show the degree of mesh stretching used in the calculations. 
It can be seen, for example, from Fig. 3 that the axial grid spacing at the aft 
edge of the slot is 0.001" while that near the outflow boundary is 10". 
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The 


radial grid distribution utilized 95 grid points while 155 grid points were 
used along the axial direction. The slot width of 1.5" was resolved using 55 
grid points while 40 grid points were used between the inflow boundary and the 
beginning of the slot, and 60 grid points were used between the aft edge of the 
slot and the outflow boundary. Figure 4 shows the computational grid used in 
the two-dimensional calculation for the entire domain in the r-z plane, while 
the details of the grid in the vicinity of the slot exit are presented in Fig. 
5. The grid for the three-dimensional calculations utilized 49 equally spaced 
points along the circumferential direction. The grid distribution in the r-z 
plane for the three-dimensional calculations was identical to the one used in 
the axisymmetric calculations. The total number of grid points for 
three-dimensional calculations was, thus, 721,525. 

Axisymmetric Flow Simulation 

The internal flow in an axisymmetric configuration was simulated to study 
the flow in the vicinity of the slot and to evaluate whether the small slot 
width along with the absence of an inhibitor on one propellant surface in the 
slot would result in the pressure at the joint being significantly above the 
mean chamber pressure. In addition, the flow field obtained for this case was 
to be used to generate the initial flow field for the three-dimensional flow 
simulation. The converged solution for the axisymmetric case was obtained in 
approximately 450 iterations, with each time step requiring approximately 2 
seconds of CPU time on a CRAY X-MP. The pressure distribution in the vicinity 
of the slot is shown in Figs. 6 and 7. Figure 6 shows the pressure 
distribution along with an outline of the computational domain while Fig. 7 
shows a close-up of the region In the vicinity of the slot. The reference 
pressure for the calculations was 900 psi and thus the over-pressure at the 
joint (i.e., location in the slot at maximum distance from the centerline) was 
only about 10 psi above the chamber pressure. The high gas velocity at the 
slot exit creates a low pressure region at the aft-edge of the slot (as can be 
seen from Fig. 8) and also creates a region of reverse axial flow near the 
wall. In Fig. 9, velocity vectors in the vicinity of the slot exit are shown 
which illustrate this phenomenon. The magnitude of the reversed axial velocity 
was found to be about 10% of the velocity at the centerline. This region of 
reversed flow was found to extend about 5 slot widths along the axial 



direction, as can be seen from Fig. 10, where streak lines in the vicinity of 
the slot exit are shown. 

Three-Dimensional Flow Simulation 


Two cases with an asymmetric flow field in the SRB chamber were simulated 
under the present contract. The objective of the three-dimensional flow 
simulation was to estimate the extent of circumferential flow and pressure 
gradients present at the bottom of the slot, as these parameters are needed in 
the analysis of the thermal/fluid environment sustained by the joint 
components. In both cases, the asymmetry in the flow was to be caused by an 
asymmetric boundary condition. In Case I, the asymmetry was assumed to be an 
asymmetric pressure imposed at the outflow boundary of the computational 
domain. This condition could exist in practice, for example, by a canted 
nozzle in the SRB. In Case II, the asymmetry was assumed to be a result of a 
failure of a part of the castable inhibitor in the slot. During the simulation 
of this case, a 7.5° section of the castable inhibitor at 0 = +90° was assumed 
to be removed, thereby generating an asymmetric flow in the slot. 

The results of the calculations for the two cases are presented in Figs. 

11 through 25 and discussed below. It should be noted that Figs. 11 through 25 
indicate the grid to be 48x49x53 rather than 95x49x155. The reason behind this 
is that, while the computations were performed on a 95x49x155 grid, the data 
had to be culled selectively to avoid exceeding the memory of the graphics 
workstation. Hence, the plot-file used consisted of approximately every 
alternate grid point in the radial direction and every third point along the 
axial direction. 

Case I Results : The asymmetric pressure distribution imposed at the outflow 

boundary was obtained by varying the static pressure at the outflow boundary 
linearly along the diameter of the outflow plane. The pressure distribution at 
the inflow boundary was axisymmetric, as can be seen from Fig. 11. The 
pressure distribution at the outflow boundary is shown in Fig. 12. The 
magnitude of the pressure drop along the diameter at the outflow plane was 
approximately 9.5 psi. The resultant asymmetry in the axial velocity can be 
seen from Fig. 13, which shows a shift in the axial velocity peak towards the 

low pressure side. The pressure and axial velocity distribution approximately 
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one diameter upstream of the outflow boundary are presented in Figs. 14 and 15, 
respectively. At this location, the axial velocity distribution appears fairly 
axisymmetric even though the pressure distribution is still not axisymmetric. 
The pressure and axial velocity distributions approximately two diameters 
upstream of the outflow boundary, shown in Figs. 16 and 17, respectively, 
indicated the flow to be axisymmetric at this location. Needless to say, the 
flow field in the slot was axisymmetric, as can be seen from the pressure 
distribution at the slot bottom shown in Fig. 18, and has no circumferential 
velocities or pressure gradients. 

Case II Results : As the results of the Case I simulation indicated no 

asymmetric effects on the flow field in the slot, the boundary condition at the 
outflow was maintained the same during Case II calculations. In addition, a 
7.5° section of the castable inhibitor at one end of the diameteral symmetry 
plane (i.e., at 0 = +90°) was assumed to be absent. The computational domain 
includes only a 180° section (-90°< 0 <+90°) of the SRB (due to symmetry), this 
corresponds to a failure of a 15° section of the castable inhibitor. The 
pressure and axial velocity at the outflow plane are shown in Figs. 19 and 20, 
respectively, which show the asymmetry imposed and the outflow boundary. The 
pressure and velocity distribution approximately two diameters upstream of the 
outflow plane, presented in Figs. 21 and 22, respectively, show that the effect 
of the asymmetry at the outflow boundary does not propagate more than two 
diameters upstream. Since, in this case, the broken castable inhibitor section 
also originated asymmetry in the slot, the flow field in the slot was observed 
to be asymmetric. The pressure distribution at the slot mid-section, shown in 
Fig. 23, indicates the maximum pressure differential along the circumferential 
direction to be approximately 1 psi. The pressure reaches a minimum value at 
0 = -30° (i.e. 120° away from the location of inhibitor failure) and a maximum 
value at a location diametrically opposite to the inhibitor failure (i.e. at 
0 = -90°). The velocity vectors in the r-0 plane at slot mid-section are shown 
in Fig. 24, with the details of the region near the inhibitor failure being 
presented in Fig. 25. 


Conclusions 


An efficient Navier-Stokes analysis has been successfully applied to 
simulate the complex flow field in the vicinity of a slot in a solid rocket 
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motor with segment joints. The capability of the computer code to resolve the 
flow near solid surfaces without using a "wall function" assumption has been 
demonstrated. In view of the complex nature of the flow field in the vicinity 
of the slot, this approach is considered essential. The results obtained from 
these calculations provide valuable design information, which would otherwise 
be extremely difficult (if not impossible) to obtain. The results of the 
axisymmetric calculations indicate the presence of a region of reversed axial 
flow at the aft-edge of the slot and show the over-pressure in the slot to be 
only about 10 psi. The results of the asymmetric calculations indicate that a 
pressure asymmetry more than two diameters downstream of the slot has no 
noticable effect on the flow field in the slot. They also indicate that the 
circumferential pressure differential caused in the slot due to failure of a 
15° section of the castable inhibitor will be approximately 1 psi. 
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Fig. 2 - Grid Distribution along Radial Direction 
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Fig. 3 - Grid Distribution along Axial Direction 
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Details of the Grid Near Slot Exit 
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Pressure Distribution in the Chamber 
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Fig. 7 - Pressure Distribution in the Vicinity of the Slot 
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Pressure Distribution Near Slot Exit 
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Velocity Vectors Near Slot Exit 
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Fig. 10 - Streaklines Near Slot Exit 
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Fig. II - Pressure Distribution at Inflow Boundary 
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Fig. 12 - Pressure Distribution at Outflow Boundary 
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Fig. 13 - Axial Velocity Distribution at Outflow Boundary 
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Fig. 14 - Pressure Distribution One Diameter Upstream of Outflow Boundary 
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Fig. 15 - Axial Velocity Distribution One Diameter Upstream of Outflow Boundary 
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Fig. 17 - Axial Velocity Distribution Two Diameters Upstream of Outflow Boundary 
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Fig. 18 - Pressure Distribution at Slot Mid-Section 
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Fig. 19 - Pressure Distribution at Outflow Boundary (Case II) 
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Fig. 20 - Axial Velocity Distribution at Outflow Boundary (Case II) 
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Fig* 21 — Pressure Distribution Two Diameters Upstream of Outflow Boundary (Case II) 
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Fig* 22 — Axial Velocity Distribution Two Diameters Upstream of Outflow Boundary (Case II) 
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Fig. 23 - Pressure Distribution at Slot Mid-Section (Case II) 
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Fig. 24 - Velocity Vectors in r-9 Plane at Slot Mid-Section (Case II 
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Fig. 25 - Velocity Vectors in r-9 Plane at Near Inhibitor Failure 



